## Replication Code for "Domestic Polarization and International Rivalry" 
## Myrick and Wang in "The Journal of Politics" 
## Last Update: February 1, 2023

## This code replicates two studies in the paper: 
## (1) survey experiment, and (2) observational analysis

# Load Packages -----------------------------------------------------------

package<- list("foreign", "MASS", "stargazer", "pastecs", "dplyr", "tidyr",
               "rmarkdown", "knitr", "ggplot2", "pscl", "ggthemes", "margins",
               "lubridate", "countrycode", "reshape", "stringr", "sandwich",
               "lmtest", "readxl", "cowplot", "latex2exp", "readstata13",
               "Zelig", "car", "writexl", "vtable", "ggpubr",
               "rstatix", "ggridges", "gridExtra", "sampleSelection")

lapply(package, FUN = function(X) {
  do.call("require", list(X)) 
})

## Note: If working in R >= 4.2, fix issue with model names in stargazer package (5.2.3) 
## using the following code

detach("package:stargazer",unload=T)
remove.packages("stargazer")
download.file("https://cran.r-project.org/src/contrib/stargazer_5.2.3.tar.gz", destfile = "stargazer_5.2.3.tar.gz")
untar("stargazer_5.2.3.tar.gz")
stargazer_src <- readLines("stargazer/R/stargazer-internal.R")
stargazer_src[1990] <- stargazer_src[1995]
stargazer_src[1995] <- ""
writeLines(stargazer_src, con="stargazer/R/stargazer-internal.R")
install.packages("stargazer", repos = NULL, type="source")
library(stargazer)

# Set Working Directory ---------------------------------------------------

## Change working directory

setwd("")

# Manuscript: Experimental Analysis ---------------------------------------------------

#Import survey data (n = 2046)

d<- read.csv("US-China Polarization - Survey Data.csv")
names(d)

#Respondents who did not pass the pre-treatment attention check
#were screened out of the survey before the experiment

d2<- filter(d, Finished=="TRUE" & attn2=="Football (Soccer)")


# Figure 1 ----------------------------------------------------------------

## Plot core test of emboldening hypothesis by looking at mean
## outcome for china_assertive DV in the control vs. treatment group

## Make china_assertive DV numeric

d2$china_assertive_num<- ifelse(d2$china_assertive=="Much less assertive", -2, NA)
d2$china_assertive_num<- ifelse(d2$china_assertive=="Somewhat less assertive", -1, d2$china_assertive_num)
d2$china_assertive_num<- ifelse(d2$china_assertive=="About the same", 0, d2$china_assertive_num)
d2$china_assertive_num<- ifelse(d2$china_assertive=="Somewhat more assertive", 1, d2$china_assertive_num)
d2$china_assertive_num<- ifelse(d2$china_assertive=="Much more assertive", 2, d2$china_assertive_num)

## Box plot + two-sample t test of china_assertive by control vs. treatment

stat.test<- t_test(china_assertive_num ~ t_polarization,
                   data=d2[which(is.na(d2$china_assertive_num)==FALSE),],
                   alternative = c("two.sided"), mu=0,
                   paired=FALSE, var.equal=FALSE, conf.level = 0.95) ## t_test is pipe-friendly and the output can be added to a plot using the ggpubr package 

d2$t_polarization<- as.factor(d2$t_polarization)

p1 <- ggplot(d2[which(is.na(d2$china_assertive_num)==FALSE),], aes(x=t_polarization, y=china_assertive_num)) + 
  geom_boxplot(width=0.5,color="grey60", fill="grey", alpha=0.2)+
  stat_summary(fun=mean, geom="point", shape=17,size=5, color="black")+
  geom_jitter(width=0.2,size=1.5, color="grey60")+
  theme_bw()+
  ylab("Preference for China's Foreign Policy")+
  xlab("")+
  scale_y_continuous(breaks = seq(-2,2,by=1), 
                     labels = c("Much less assertive (-2)", 
                                "Somewhat less assertive (-1)",
                                "About the same (0)", 
                                "Somewhat more assertive (1)", 
                                "Much more assertive (2)"))+
  scale_x_discrete()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 15, b = 0, l = 0), size=20),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0), size=20),
        text = element_text(size=25),
        legend.text=element_text(size=25),
        legend.key.size=unit(0.8, "cm"),
        legend.position="")+
  labs(subtitle = get_test_label(stat.test, detailed = TRUE))

p1

# Figure 2 ----------------------------------------------------------------

## Plot relative importance of foreign policy issues to China vs. US
## to identify issues core to China's national interests

d3<- select(d2, ResponseId,china_importance_1, china_importance_2, china_importance_3,
            china_importance_4, china_importance_5, china_importance_6, 
            us_importance_1, us_importance_2, us_importance_3,
            us_importance_4, us_importance_5, us_importance_6,
            t_polarization)

d3<- dplyr::rename(d3, tw_c=china_importance_1, scs_c=china_importance_2,
                   cs_c=china_importance_3, tsc_c=china_importance_4,
                   gl_c=china_importance_5, os_c=china_importance_6,
                   tw_u=us_importance_1, scs_u=us_importance_2,
                   cs_u=us_importance_3, tsc_u=us_importance_4,
                   gl_u=us_importance_5, os_u=us_importance_6)

## As specified in PAP, use only data from control group 
dd_control<- d3[which(d3$t_polarization=="Control"),]

## Paired t-tests for relative importance (China vs. US)
issue_lista<- list("tw_c", "scs_c", "cs_c", "tsc_c", "gl_c", "os_c")

issue_listb<- list("tw_u", "scs_u", "cs_u", "tsc_u", "gl_u", "os_u")


result<- data.frame(matrix(nrow=6, ncol=3))
colnames(result) <- c("meandiff", "lower", "upper")
number<- c(1, 2, 3, 4, 5, 6)


for (j in number) {
  dtest<- select(dd_control, ResponseId, issue_lista[[j]], issue_listb[[j]])
  t<- t.test(dtest[,2], dtest[,3], alternative = c("two.sided"), mu=0,
              paired=TRUE, conf.level = 0.95)
  result[j, 1]<- t$estimate[[1]]
  result[j, 2]<- t$conf.int[[1]]
  result[j, 3]<- t$conf.int[[2]]
}

result$type<- as.factor(c(1,2,3,4,5,6))
result$cat<- c("Taiwan", "South China Sea", "Cybersecurity", 
               "Trade & Supply Chain", "Global Leadership", "Outerspace")



p2 <- ggplot(result, aes(x=type, y=meandiff)) + 
  geom_point(size=5)+
  coord_flip()+
  geom_errorbar(aes(ymin = lower, ymax=upper),size=1.2, width=0.1) +
  geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  theme_bw()+
  ylab("Relative importance")+
  xlab("")+
  ggtitle("(a).Relative Issue Salience (China vs. U.S.)")+
  scale_x_discrete(labels = c("Taiwan", "South \n China Sea", "Cybersecurity", 
                              "Trade & \n Supply Chain", "Global Leadership", "Outerspace"))+
  theme(axis.title.y = element_text(margin = margin(t = 20, r = 15, b = 0, l = 0), size=30, face="bold"),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0), size=30, face="bold"),
        text = element_text(size=35, face="bold"),
        plot.title = element_text(face = "bold", size=35),
        legend.text=element_text(size=25),
        legend.key.size=unit(0.8, "cm"),
        legend.position="top")

p2

## Plot t tests for issue-specific assertiveness

#### Make six DVs (issue-specific) numeric
#### On Taiwan
d2$china_issue_1_num<- ifelse(d2$china_issue_1=="Much less assertive", -2, NA)
d2$china_issue_1_num<- ifelse(d2$china_issue_1=="Somewhat less assertive", -1, d2$china_issue_1_num)
d2$china_issue_1_num<- ifelse(d2$china_issue_1=="About the Same", 0, d2$china_issue_1_num)
d2$china_issue_1_num<- ifelse(d2$china_issue_1=="Somewhat more assertive", 1, d2$china_issue_1_num)
d2$china_issue_1_num<- ifelse(d2$china_issue_1=="Much more assertive", 2, d2$china_issue_1_num)

#### On South China Sea
d2$china_issue_2_num<- ifelse(d2$china_issue_2=="Much less assertive", -2, NA)
d2$china_issue_2_num<- ifelse(d2$china_issue_2=="Somewhat less assertive", -1, d2$china_issue_2_num)
d2$china_issue_2_num<- ifelse(d2$china_issue_2=="About the Same", 0, d2$china_issue_2_num)
d2$china_issue_2_num<- ifelse(d2$china_issue_2=="Somewhat more assertive", 1, d2$china_issue_2_num)
d2$china_issue_2_num<- ifelse(d2$china_issue_2=="Much more assertive", 2, d2$china_issue_2_num)

#### On Cybersecurity
d2$china_issue_3_num<- ifelse(d2$china_issue_3=="Much less assertive", -2, NA)
d2$china_issue_3_num<- ifelse(d2$china_issue_3=="Somewhat less assertive", -1, d2$china_issue_3_num)
d2$china_issue_3_num<- ifelse(d2$china_issue_3=="About the Same", 0, d2$china_issue_3_num)
d2$china_issue_3_num<- ifelse(d2$china_issue_3=="Somewhat more assertive", 1, d2$china_issue_3_num)
d2$china_issue_3_num<- ifelse(d2$china_issue_3=="Much more assertive", 2, d2$china_issue_3_num)

#### On Trade and Supply Chain
d2$china_issue_4_num<- ifelse(d2$china_issue_4=="Much less assertive", -2, NA)
d2$china_issue_4_num<- ifelse(d2$china_issue_4=="Somewhat less assertive", -1, d2$china_issue_4_num)
d2$china_issue_4_num<- ifelse(d2$china_issue_4=="About the Same", 0, d2$china_issue_4_num)
d2$china_issue_4_num<- ifelse(d2$china_issue_4=="Somewhat more assertive", 1, d2$china_issue_4_num)
d2$china_issue_4_num<- ifelse(d2$china_issue_4=="Much more assertive", 2, d2$china_issue_4_num)


#### On Global Leadership
d2$china_issue_5_num<- ifelse(d2$china_issue_5=="Much less assertive", -2, NA)
d2$china_issue_5_num<- ifelse(d2$china_issue_5=="Somewhat less assertive", -1, d2$china_issue_5_num)
d2$china_issue_5_num<- ifelse(d2$china_issue_5=="About the Same", 0, d2$china_issue_5_num)
d2$china_issue_5_num<- ifelse(d2$china_issue_5=="Somewhat more assertive", 1, d2$china_issue_5_num)
d2$china_issue_5_num<- ifelse(d2$china_issue_5=="Much more assertive", 2, d2$china_issue_5_num)

#### On Outerspace
d2$china_issue_6_num<- ifelse(d2$china_issue_6=="Much less assertive", -2, NA)
d2$china_issue_6_num<- ifelse(d2$china_issue_6=="Somewhat less assertive", -1, d2$china_issue_6_num)
d2$china_issue_6_num<- ifelse(d2$china_issue_6=="About the Same", 0, d2$china_issue_6_num)
d2$china_issue_6_num<- ifelse(d2$china_issue_6=="Somewhat more assertive", 1, d2$china_issue_6_num)
d2$china_issue_6_num<- ifelse(d2$china_issue_6=="Much more assertive", 2, d2$china_issue_6_num)

issue_response_list<- list("china_issue_1_num", "china_issue_2_num", 
                           "china_issue_3_num", "china_issue_4_num", 
                           "china_issue_5_num", "china_issue_6_num")

result2<- data.frame(matrix(nrow=6, ncol=3))
colnames(result2) <- c("meandiff", "lower", "upper")
number<- c(1, 2, 3, 4, 5, 6)

for (j in number) {
  dtest<- select(d2, t_polarization, issue_response_list[[j]])
  names(dtest)[2]<- "dv"
  t<- t.test(dtest[which(dtest$t_polarization=="Treatment"),]$dv,
             dtest[which(dtest$t_polarization=="Control" ),]$dv,
             alternative = c("two.sided"), mu=0,
             paired=FALSE, var.equal=FALSE, conf.level = 0.95)
  result2[j, 1]<- t$estimate[[1]]-t$estimate[[2]]
  result2[j, 2]<- t$conf.int[[1]]
  result2[j, 3]<- t$conf.int[[2]]
}

result2$type<- as.factor(c(1,2,3,4,5,6))
result2$cat<- c("Taiwan", "South China Sea", "Cybersecurity", 
                "Trade & Supply Chain", "Global Leadership", "Outerspace")


p3 <- ggplot(result2, aes(x=type, y=meandiff)) + 
  geom_point(size=5)+
  coord_flip()+
  geom_errorbar(aes(ymin = lower, ymax=upper),size=1.2, width=0.1) +
  geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  theme_bw()+
  ylab("Difference in Means \n (Treatment - Control)")+
  xlab("")+
  ylim(-0.5, 0.5)+
  ggtitle("(b).Selective Assertiveness (Treatment vs. Control)")+
  scale_x_discrete(labels = c("Taiwan", "South \n China Sea", "Cybersecurity", 
                              "Trade & \n Supply Chain", "Global Leadership", "Outerspace"))+
  theme(axis.title.y = element_text(margin = margin(t = 20, r = 15, b = 0, l = 0), size=30, face="bold"),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0), size=30, face="bold"),
        text = element_text(size=35, face="bold"),
        plot.title = element_text(face = "bold", size=35),
        legend.text=element_text(size=25),
        legend.key.size=unit(0.8, "cm"),
        legend.position="top")

p3

# Figure 2

cowplot::plot_grid(p2, p3, labels = "", ncol = 2, 
                   align = 'v', axis = 'l') 


# Figure 3 ----------------------------------------------------------------

## Plot mean outcomes for US foreign policy (us_strength, us_active, us_assertive)
## by control vs. treatment group

#### Make US DVs numeric

d2$us_strength_num<- ifelse(d2$us_strength=="Much weaker", -2, NA)
d2$us_strength_num<- ifelse(d2$us_strength=="Somewhat weaker", -1, d2$us_strength_num)
d2$us_strength_num<- ifelse(d2$us_strength=="Neither weaker nor stronger", 0, d2$us_strength_num)
d2$us_strength_num<- ifelse(d2$us_strength=="Somewhat stronger", 1, d2$us_strength_num)
d2$us_strength_num<- ifelse(d2$us_strength=="Much stronger", 2, d2$us_strength_num)

d2$us_active_num<- ifelse(d2$us_active=="Much less active in global affairs", -2, NA)
d2$us_active_num<- ifelse(d2$us_active=="Somewhat less active in global affairs", -1, d2$us_active_num)
d2$us_active_num<- ifelse(d2$us_active=="About the same", 0, d2$us_active_num)
d2$us_active_num<- ifelse(d2$us_active=="Somewhat more active in global affairs", 1, d2$us_active_num)
d2$us_active_num<- ifelse(d2$us_active=="Much more active in global affairs", 2, d2$us_active_num)

d2$us_assertive_num<- ifelse(d2$us_assertive=="Much less assertive", -2, NA)
d2$us_assertive_num<- ifelse(d2$us_assertive=="Somewhat less assertive", -1, d2$us_assertive_num)
d2$us_assertive_num<- ifelse(d2$us_assertive=="About the same", 0, d2$us_assertive_num)
d2$us_assertive_num<- ifelse(d2$us_assertive=="Somewhat more assertive", 1, d2$us_assertive_num)
d2$us_assertive_num<- ifelse(d2$us_assertive=="Much more assertive", 2, d2$us_assertive_num)

## US_STRENGTH

t.test(us_strength_num ~ t_polarization,
       data=d2[which(is.na(d2$us_strength_num)==FALSE),],
       alternative = c("two.sided"), mu=0,
       paired=FALSE, var.equal=FALSE, conf.level = 0.95)

stat.test<- t_test(us_strength_num ~ t_polarization,
                   data=d2[which(is.na(d2$us_strength_num)==FALSE),],
                   alternative = c("two.sided"), mu=0,
                   paired=FALSE, var.equal=FALSE, conf.level = 0.95) ## t_test is pipe-friendly and the output can be added to a plot using the ggpubr package 

d2$t_polarization<- as.factor(d2$t_polarization)

p4 <- ggplot(d2[which(is.na(d2$us_strength_num)==FALSE),], aes(x=t_polarization, y=us_strength_num)) + 
  geom_boxplot(width=0.5,color="grey60", fill="grey", alpha=0.2)+
  stat_summary(fun=mean, geom="point", shape=17,size=5, color="black")+
  geom_jitter(width=0.2,size=0.6, color="grey60")+
  theme_bw()+
  ylab("")+
  xlab("")+
  ggtitle("(a) Perception of US strength")+
  scale_y_continuous(breaks = seq(-2,2,by=1), 
                     labels = c("Much weaker (-2)", 
                                "Somewhat weaker (-1)",
                                "Neither weaker \n nor stronger (0)", 
                                "Somewhat stronger (1)", 
                                "Much stronger (2)"))+
  scale_x_discrete()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 5, b = 0, l = 0), size=20),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0), size=20),
        text = element_text(size=35, face="bold"),
        plot.title = element_text(face = "bold", size=30),
        plot.subtitle = element_text(vjust = -0.5),
        legend.text=element_text(size=25),
        legend.key.size=unit(0.8, "cm"),
        legend.position="",
        plot.margin = unit(c(0.4, 0.4, 0.4, 0.4), "cm"))+
  labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p4

## US_ACTIVE

t.test(us_active_num ~ t_polarization,
       data=d2[which(is.na(d2$us_active_num)==FALSE),],
       alternative = c("two.sided"), mu=0,
       paired=FALSE, var.equal=FALSE, conf.level = 0.95) 

stat.test<- t_test(us_active_num ~ t_polarization,
                   data=d2[which(is.na(d2$us_active_num)==FALSE),],
                   alternative = c("two.sided"), mu=0,
                   paired=FALSE, var.equal=FALSE, conf.level = 0.95) ## t_test is pipe-friendly and the output can be added to a plot using the ggpubr package 

d2$t_polarization<- as.factor(d2$t_polarization)

p5 <- ggplot(d2[which(is.na(d2$us_active_num)==FALSE),], aes(x=t_polarization, y=us_active_num)) + 
  geom_boxplot(width=0.5,color="grey60", fill="grey", alpha=0.2)+
  stat_summary(fun=mean, geom="point", shape=17,size=5, color="black")+
  geom_jitter(width=0.2,size=0.6, color="grey60")+
  theme_bw()+
  ylab("")+
  xlab("")+
  ggtitle("(b) Perception of US activeness in global affairs")+
  scale_y_continuous(breaks = seq(-2,2,by=1), 
                     labels = c("Much less active (-2)", 
                                "Somewhat less active (-1)",
                                "About the same (0)", 
                                "Somewhat more active (1)", 
                                "Much more active (2)"))+
  scale_x_discrete()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 5, b = 0, l = 0), size=20),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0), size=20),
        text = element_text(size=35, face="bold"),
        plot.title = element_text(face = "bold", size=30),
        plot.subtitle = element_text(vjust = -0.5),
        legend.text=element_text(size=30),
        legend.key.size=unit(0.8, "cm"),
        legend.position="",
        plot.margin = unit(c(0.4, 0.4, 0.4, 0.4), "cm"))+
  labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p5


## US_ASSERTIVE

t.test(us_assertive_num ~ t_polarization,
       data=d2[which(is.na(d2$us_assertive_num)==FALSE),],
       alternative = c("two.sided"), mu=0,
       paired=FALSE, var.equal=FALSE, conf.level = 0.95)

stat.test<- t_test(us_assertive_num ~ t_polarization,
                   data=d2[which(is.na(d2$us_assertive_num)==FALSE),],
                   alternative = c("two.sided"), mu=0,
                   paired=FALSE, var.equal=FALSE, conf.level = 0.95) ## t_test is pipe-friendly and the output can be added to a plot using the ggpubr package 

d2$t_polarization<- as.factor(d2$t_polarization)

p6 <- ggplot(d2[which(is.na(d2$us_assertive_num)==FALSE),], aes(x=t_polarization, y=us_assertive_num)) + 
  geom_boxplot(width=0.5,color="grey60", fill="grey", alpha=0.2)+
  stat_summary(fun=mean, geom="point", shape=17,size=5, color="black")+
  geom_jitter(width=0.2,size=0.6, color="grey60")+
  theme_bw()+
  ylab("")+
  xlab("")+
  ggtitle("(c) Perception of US assertiveness towards China")+
  scale_y_continuous(breaks = seq(-2,2,by=1), 
                     labels = c("Much less assertive (-2)", 
                                "Somewhat less assertive (-1)",
                                "About the same (0)", 
                                "Somewhat more assertive (1)", 
                                "Much more assertive (2)"))+
  scale_x_discrete(labels=c("Control Group", "Treatment Group"))+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 5, b = 0, l = 0), size=20),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0), size=20),
        text = element_text(size=35, face="bold"),
        plot.title = element_text(face = "bold", size=30),
        plot.subtitle = element_text(vjust = -0.5),
        legend.text=element_text(size=30),
        legend.key.size=unit(0.8, "cm"),
        legend.position="",
        plot.margin = unit(c(0.4, 0.4, 0.4, 0.4), "cm"))+
  labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p6


# Figure 3

library(gridExtra)
grid.arrange(p4, p5, p6, 
             ncol = 2, 
             nrow = 2,
             layout_matrix = cbind(c(1,1), c(2,3)))


# Appendix: Experimental Analysis -----------------------------------------


# Appendix I-A: Demographics ------------------------------------------------------------
rm(list = ls())

#Import survey data (n = 2046)

d<- read.csv("US-China Polarization - Survey Data.csv")
names(d)

#Respondents who did not pass the pre-treatment attention check
#were screened out of the survey before the experiment

d2<- filter(d, Finished=="TRUE" & attn2=="Football (Soccer)")

## Code sex of respondent
d2$female = ifelse(d2$sex=="Female", 1, 0)
d2$female = ifelse(d2$sex=="", NA, d2$female)

## Code age of respondent
d2$age = as.numeric(d2$age)

# Bin age into three dummy variables as specified in PAP
# Code incorrect age entries as missing
d2$age1834 = ifelse(d2$age > 100 | d2$age < 18, NA, d2$age) 
d2$age1834 = ifelse(d2$age1834 > 17 & d2$age1834 < 35, 1, 0)
d2$age1834 = ifelse(is.na(d2$age)==TRUE, NA, d2$age1834)
d2$age3554 = ifelse(d2$age > 100 | d2$age < 18, NA, d2$age)
d2$age3554 = ifelse(d2$age3554 > 34 & d2$age3554 < 55, 1, 0)
d2$age3554 = ifelse(is.na(d2$age)==TRUE, NA, d2$age3554)
d2$age55plus = ifelse(d2$age > 100 | d2$age < 18, NA, d2$age)
d2$age55plus = ifelse(d2$age55plus > 54, 1, 0)
d2$age55plus = ifelse(is.na(d2$age)==TRUE, NA, d2$age55plus)

# Code Han ethnicity
d2$han = ifelse(d2$ethnicity=="Han", 1, 0)
d2$han = ifelse(d2$ethnicity=="", NA, d2$han)

# Code Urban/Rural 
d2$urban_bi = ifelse(d2$urban == "Urban", 1, 0)
d2$urban_bi = ifelse(d2$urban == "", NA, d2$urban_bi)

# Code Education
d2$college_higher = ifelse(d2$education == "College (including Junior College)" |
                             d2$education == "Doctoral" |
                             d2$education == "Master", 1, 0)
d2$college_higher = ifelse(d2$education == "", NA, d2$college_higher)

## Compare Survey Resposes to Target Demographic Quotas

  # Age (Target: 30% 18-34 / 32% 35-54 / 38% 55+)
prop.table(table(d2$age1834))
prop.table(table(d2$age3554))
prop.table(table(d2$age55plus))

  # Sex (Target: 48% M / 52% F)
prop.table(table(d2$female))

## Look at other Demographic Information

  # Urban Population
prop.table(table(d2$urban_bi)) #97.3%
  
  # Ethnicity (Percent Han)
prop.table(table(d2$han)) #96.7%

  # Education
prop.table(table(d2$college_higher)) #78.6%

 # region
prop.table(table(d2$region))


# Appendix I-C: Sample Quality --------------------------------------------------------

d2<- filter(d, Finished=="TRUE" & attn2=="Football (Soccer)")

## Time spent reading polarization prime: Figure 1 in Appendix

dtreatment<- filter(d2, t_polarization=="Treatment")
dtreatment<- select(dtreatment, ResponseId, Q54_Page.Submit, Q55_Page.Submit,
                    Q56_Page.Submit, Q57_Page.Submit)
dtreatment$prime_time = dtreatment$Q54_Page.Submit + dtreatment$Q55_Page.Submit +
  dtreatment$Q56_Page.Submit + dtreatment$Q57_Page.Submit

prime_time <- ggplot(dtreatment, aes(x=prime_time)) + theme_bw(base_size = 25)+
  geom_density(fill="gray") +  ylab("Density")+
  xlab("Time (in Seconds)")+   ggtitle("Time Spent Reading Prime in Treatment Group") +
  geom_vline(xintercept = median(dtreatment$prime_time, na.rm = TRUE), color = "red", linetype="dashed") + 
  scale_x_continuous(limits=c(0, 150)) 

prime_time

## Manipulation Check: Table 1 in Appendix

d2$manipulation_check_bi1<- ifelse(d2$manipulation_check=="Almost Always" | d2$manipulation_check=="Sometimes", 1, 0)
d2$manipulation_check_bi1<- ifelse(d2$manipulation_check=="I don't know", NA, d2$manipulation_check_bi1)
d2$manipulation_check_bi1<- ifelse(d2$manipulation_check=="", NA, d2$manipulation_check_bi1) 


m_check = glm(manipulation_check_bi1 ~ t_polarization, data =d2, family = binomial(link = "logit"))

stargazer(m_check, title = "How Often Would You Say These Parties Agree?",
          dep.var.labels = "Almost Always/Sometimes",
          covariate.labels = "Treatment")

# Appendix I-D: Interpretation of "Assertiveness" ----------------------------------------------------------

d2$expand_defend<- ifelse(str_detect(d2$assertive_def, "Defend"), 0, NA)
d2$expand_defend<- ifelse(str_detect(d2$assertive_def, "Expand"), 1, d2$expand_defend)

d2$res_act_agg<- ifelse(str_detect(d2$assertive_def, "resolutely"), 0, NA)
d2$res_act_agg<- ifelse(str_detect(d2$assertive_def, "actively"), 1, d2$res_act_agg)
d2$res_act_agg<- ifelse(str_detect(d2$assertive_def, "aggressively"), 2, d2$res_act_agg)

p_interpretation <- ggplot(d2, aes(x=expand_defend, y=res_act_agg)) + 
  geom_point(width=0.5)+
  #stat_summary(fun=mean, geom="point", size=4, color="red")+
  geom_jitter(width=0.3,size=2)+
  theme_bw()+
  ylab("")+
  xlab("")+
  scale_y_continuous(breaks = seq(0,2,by=1), 
                     labels = c("Resolutely", 
                                "Actively",
                                "Aggressively"))+
  scale_x_continuous(breaks = seq(0,1, by=1),
                     labels = c("Defend China's interests",
                                "Expand China's influence"))+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 15, b = 0, l = 0), size=20),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0), size=20),
        text = element_text(size=25, face="bold"),
        legend.text=element_text(size=25),
        legend.key.size=unit(0.8, "cm"),
        legend.position="")

# Figure 2 in Appendix

p_interpretation

# Appendix I-E: Robustness Checks for Survey Experiment ------------------------------

# Table 2 in Appendix: Regression Analysis of Preference for Chinese Assertiveness
d2$polarization_treatment = ifelse(d2$t_polarization=="Treatment", 1, 0)

## Make DV numeric
d2$china_assertive_num<- ifelse(d2$china_assertive=="Much less assertive", -2, NA)
d2$china_assertive_num<- ifelse(d2$china_assertive=="Somewhat less assertive", -1, d2$china_assertive_num)
d2$china_assertive_num<- ifelse(d2$china_assertive=="About the same", 0, d2$china_assertive_num)
d2$china_assertive_num<- ifelse(d2$china_assertive=="Somewhat more assertive", 1, d2$china_assertive_num)
d2$china_assertive_num<- ifelse(d2$china_assertive=="Much more assertive", 2, d2$china_assertive_num)

## Code demographic controls
d2$female = ifelse(d2$sex=="Female", 1, 0)
d2$female = ifelse(d2$sex=="", NA, d2$female)

d2$age = as.numeric(d2$age)

d2$age1834 = ifelse(d2$age > 100 | d2$age < 18, NA, d2$age)
d2$age1834 = ifelse(d2$age1834 > 17 & d2$age1834 < 35, 1, 0)
d2$age1834 = ifelse(is.na(d2$age)==TRUE, NA, d2$age1834)
d2$age3554 = ifelse(d2$age > 100 | d2$age < 18, NA, d2$age)
d2$age3554 = ifelse(d2$age3554 > 34 & d2$age3554 < 55, 1, 0)
d2$age3554 = ifelse(is.na(d2$age)==TRUE, NA, d2$age3554)
d2$age55plus = ifelse(d2$age > 100 | d2$age < 18, NA, d2$age)
d2$age55plus = ifelse(d2$age55plus > 54, 1, 0)
d2$age55plus = ifelse(is.na(d2$age)==TRUE, NA, d2$age55plus)

d2$han = ifelse(d2$ethnicity=="Han", 1, 0)
d2$han = ifelse(d2$ethnicity=="", NA, d2$han)
d2$urban_bi = ifelse(d2$urban == "Urban", 1, 0)
d2$urban_bi = ifelse(d2$urban == "", NA, d2$urban_bi)
d2$college_higher = ifelse(d2$education == "College (including Junior College)" |
                             d2$education == "Doctoral" |
                             d2$education == "Master", 1, 0)
d2$college_higher = ifelse(d2$education == "", NA, d2$college_higher)

m1<- lm(china_assertive_num ~ polarization_treatment, data = d2)
summary(m1)

m2<- lm(china_assertive_num ~ polarization_treatment  + female + age3554 + 
          age55plus + han + urban_bi + college_higher, data = d2)
summary(m2)

d2$china_assertive_bi<- ifelse(d2$china_assertive_num>=1, 1, 0)
d2$china_assertive_bi<- ifelse(d2$china_assertive_num=="", NA, d2$china_assertive_bi)

m3<- glm(china_assertive_bi~ polarization_treatment, family = "binomial",
         data=d2)
summary(m3)

m4<- glm(china_assertive_bi~ polarization_treatment  + female + age3554 + 
           age55plus + han + urban_bi + college_higher, family = "binomial",
         data=d2)
summary(m4)

stargazer(m1, m2, m3, m4,
          align=TRUE, no.space = TRUE,
          column.sep.width = "-16pt",
          covariate.labels=c("Treatment", "Female", "Age (35-54)",
                             "Age (55 Plus)", "Han Ethnic", "Urban",
                             "Higher Education"),
          dep.var.caption = "Assertiveness of China",
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001")

## Table 3: Regression Analysis of Perceptions of the U.S. Strength and Foreign Policy

d2$us_strength_num<- ifelse(d2$us_strength=="Much weaker", -2, NA)
d2$us_strength_num<- ifelse(d2$us_strength=="Somewhat weaker", -1, d2$us_strength_num)
d2$us_strength_num<- ifelse(d2$us_strength=="Neither weaker nor stronger", 0, d2$us_strength_num)
d2$us_strength_num<- ifelse(d2$us_strength=="Somewhat stronger", 1, d2$us_strength_num)
d2$us_strength_num<- ifelse(d2$us_strength=="Much stronger", 2, d2$us_strength_num)

d2$us_active_num<- ifelse(d2$us_active=="Much less active in global affairs", -2, NA)
d2$us_active_num<- ifelse(d2$us_active=="Somewhat less active in global affairs", -1, d2$us_active_num)
d2$us_active_num<- ifelse(d2$us_active=="About the same", 0, d2$us_active_num)
d2$us_active_num<- ifelse(d2$us_active=="Somewhat more active in global affairs", 1, d2$us_active_num)
d2$us_active_num<- ifelse(d2$us_active=="Much more active in global affairs", 2, d2$us_active_num)

d2$us_assertive_num<- ifelse(d2$us_assertive=="Much less assertive", -2, NA)
d2$us_assertive_num<- ifelse(d2$us_assertive=="Somewhat less assertive", -1, d2$us_assertive_num)
d2$us_assertive_num<- ifelse(d2$us_assertive=="About the same", 0, d2$us_assertive_num)
d2$us_assertive_num<- ifelse(d2$us_assertive=="Somewhat more assertive", 1, d2$us_assertive_num)
d2$us_assertive_num<- ifelse(d2$us_assertive=="Much more assertive", 2, d2$us_assertive_num)

m1<- lm(us_strength_num ~ polarization_treatment, data = d2)
summary(m1)

m2<- lm(us_strength_num ~ polarization_treatment  + female + age3554 + 
          age55plus + han + urban_bi + college_higher, data = d2)
summary(m2)

m3<- lm(us_active_num ~ polarization_treatment, data = d2)
summary(m3)

m4<- lm(us_active_num ~ polarization_treatment  + female + age3554 + 
          age55plus + han + urban_bi + college_higher, data = d2)
summary(m4)

m5<- lm(us_assertive_num ~ polarization_treatment, data = d2)
summary(m5)

m6<- lm(us_assertive_num ~ polarization_treatment  + female + age3554 + 
          age55plus + han + urban_bi + college_higher, data = d2)
summary(m6)


stargazer(m1, m2, m3, m4, m5, m6,
          align=TRUE, no.space = TRUE,
          column.sep.width = "-16pt",
          covariate.labels=c("Treatment", "Female", "Age (35-54)",
                             "Age (55 Plus)", "Han Ethnic", "Urban",
                             "Higher Education"),
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001"
)

# Appendix I-F: Relative Importance Plots ------------------------------

d3<- select(d2, ResponseId,china_importance_1, china_importance_2, china_importance_3,
            china_importance_4, china_importance_5, china_importance_6, 
            us_importance_1, us_importance_2, us_importance_3,
            us_importance_4, us_importance_5, us_importance_6,
            t_polarization)

d3<- dplyr::rename(d3, tw_c=china_importance_1, scs_c=china_importance_2,
                   cs_c=china_importance_3, tsc_c=china_importance_4,
                   gl_c=china_importance_5, os_c=china_importance_6,
                   tw_u=us_importance_1, scs_u=us_importance_2,
                   cs_u=us_importance_3, tsc_u=us_importance_4,
                   gl_u=us_importance_5, os_u=us_importance_6)

d3$tw<- d3$tw_c-d3$tw_u
d3$scs<- d3$scs_c-d3$scs_u
d3$cs<- d3$cs_c-d3$cs_u
d3$tsc<- d3$tsc_c-d3$tsc_u
d3$gl<- d3$gl_c-d3$gl_u
d3$os<- d3$os_c-d3$os_u

dd<- select(d3, ResponseId, tw, scs, cs, tsc, gl, os, t_polarization)
dd<- gather(dd, issue, relative_importance, 2:7, factor_key = TRUE) 

dd_control<- dd[which(dd$t_polarization=="Control"),]
dd_treatment<- dd[which(dd$t_polarization=="Treatment"),]

issue_list<- list("tw", "scs", "cs", "tsc", "gl", "os")

# Relative Importance Plot: Control Group

result<- data.frame(matrix(nrow=6, ncol=3))
colnames(result) <- c("meandiff", "lower", "upper")
number<- c(1, 2, 3, 4, 5, 6)

for (j in number) {
  dtest<- filter(dd_control, issue==issue_list[[j]])
  t<- t.test(dtest$relative_importance, alternative = c("two.sided"), mu=0,
             paired=FALSE, var.equal=FALSE, conf.level = 0.95)
  result[j, 1]<- t$estimate[[1]]
  result[j, 2]<- t$conf.int[[1]]
  result[j, 3]<- t$conf.int[[2]]
}

result$type<- as.factor(c(1,2,3,4,5,6))
result$cat<- c("Taiwan", "South China Sea", "Cybersecurity", 
               "Trade & Supply Chain", "Global Leadership", "Outerspace")


p2_control <- ggplot(result, aes(x=type, y=meandiff)) + 
  geom_point(size=5)+
  geom_errorbar(aes(ymin = lower, ymax=upper),size=1.2, width=0.1) +
  geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  theme_bw()+
  ylab("")+
  xlab("")+
  ggtitle("(a). Control Group")+
  scale_x_discrete(labels = c("Taiwan", "South \n China Sea", "Cybersecurity", 
                              "Trade & \n Supply Chain", "Global Leadership", "Outerspace"))+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 15, b = 0, l = 0), size=20),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0), size=20),
        text = element_text(size=25),
        plot.title = element_text(face = "bold", size=25),
        legend.text=element_text(size=25),
        legend.key.size=unit(0.8, "cm"),
        legend.position="top")

p2_control

# Relative Importance Plot: Treatment Group

result<- data.frame(matrix(nrow=6, ncol=3))
colnames(result) <- c("meandiff", "lower", "upper")
number<- c(1, 2, 3, 4, 5, 6)

for (j in number) {
  dtest<- filter(dd_treatment, issue==issue_list[[j]])
  t<- t.test(dtest$relative_importance, alternative = c("two.sided"), mu=0,
             paired=FALSE, var.equal=FALSE, conf.level = 0.95)
  result[j, 1]<- t$estimate[[1]]
  result[j, 2]<- t$conf.int[[1]]
  result[j, 3]<- t$conf.int[[2]]
}

result$type<- as.factor(c(1,2,3,4,5,6))
result$cat<- c("Taiwan", "South China Sea", "Cybersecurity", 
               "Trade & Supply Chain", "Global Leadership", "Outerspace")


p2_treatment <- ggplot(result, aes(x=type, y=meandiff)) + 
  geom_point(size=5)+
  geom_errorbar(aes(ymin = lower, ymax=upper),size=1.2, width=0.1) +
  geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  theme_bw()+
  ylab("")+
  xlab("")+
  ggtitle("(b). Treatment Group")+
  scale_x_discrete(labels = c("Taiwan", "South \n China Sea", "Cybersecurity", 
                              "Trade & \n Supply Chain", "Global Leadership", "Outerspace"))+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 15, b = 0, l = 0), size=20),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0), size=20),
        text = element_text(size=25),
        plot.title = element_text(face = "bold", size=25),
        legend.text=element_text(size=25),
        legend.key.size=unit(0.8, "cm"),
        legend.position="top")

p2_treatment

# Relative Importance Plot: Full Sample

result<- data.frame(matrix(nrow=6, ncol=3))
colnames(result) <- c("meandiff", "lower", "upper")
number<- c(1, 2, 3, 4, 5, 6)

for (j in number) {
  dtest<- filter(dd, issue==issue_list[[j]])
  t<- t.test(dtest$relative_importance, alternative = c("two.sided"), mu=0,
             paired=FALSE, var.equal=FALSE, conf.level = 0.95)
  result[j, 1]<- t$estimate[[1]]
  result[j, 2]<- t$conf.int[[1]]
  result[j, 3]<- t$conf.int[[2]]
}

result$type<- as.factor(c(1,2,3,4,5,6))
result$cat<- c("Taiwan", "South China Sea", "Cybersecurity", 
               "Trade & Supply Chain", "Global Leadership", "Outerspace")


p2_full <- ggplot(result, aes(x=type, y=meandiff)) + 
  geom_point(size=5)+
  geom_errorbar(aes(ymin = lower, ymax=upper),size=1.2, width=0.1) +
  geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  theme_bw()+
  ylab("")+
  xlab("")+
  ggtitle("(c). Full Sample")+
  scale_x_discrete(labels = c("Taiwan", "South \n China Sea", "Cybersecurity", 
                              "Trade & \n Supply Chain", "Global Leadership", "Outerspace"))+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 15, b = 0, l = 0), size=20),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0), size=20),
        text = element_text(size=25),
        plot.title = element_text(face = "bold", size=25),
        legend.text=element_text(size=25),
        legend.key.size=unit(0.8, "cm"),
        legend.position="top")

p2_full

# Appendix Figure 3

pc<- grid.arrange(p2_control, p2_treatment, p2_full, 
                  ncol = 1, 
                  nrow = 3)
annotate_figure(pc, left = text_grob("Relative importance (one sample t-test)", 
                                     rot = 90, vjust = 0.5, 
                                     face = "bold", size = 30))

# Appendix I-G: Perceptions of Volatility in U.S. Foreign Policy --------------------------

##US_ACTIVE_CONF
stat.test<- t_test(us_active_conf_1 ~ t_polarization,
                   data=d2[which(is.na(d2$us_active_conf_1)==FALSE),],
                   alternative = c("two.sided"), mu=0,
                   paired=FALSE, var.equal=FALSE, conf.level = 0.95) ## t_test is pipe-friendly and the output can be added to a plot using the ggpubr package 

d2$t_polarization<- as.factor(d2$t_polarization)

p1 <- ggplot(d2[which(is.na(d2$us_active_conf_1)==FALSE),], aes(x=t_polarization, y=us_active_conf_1)) + 
  geom_boxplot(width=0.5,color="grey60", fill="grey", alpha=0.2)+
  stat_summary(fun=mean, geom="point", shape=17,size=5, color="black")+
  geom_jitter(width=0.2,size=1, color="grey60")+
  theme_bw()+
  ylab("Confidence about predictations \n of US activeness")+
  xlab("")+
  scale_y_continuous(breaks = seq(0,100,by=10))+
  scale_x_discrete()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 15, b = 0, l = 0), size=30, face="bold"),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0), size=20),
        text = element_text(size=30,face="bold"),
        legend.text=element_text(size=25),
        legend.key.size=unit(0.8, "cm"),
        legend.position="")+
  labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p1

##US_ASSERTIVE_CONF
stat.test<- t_test(us_assertive_conf_1 ~ t_polarization,
                   data=d2[which(is.na(d2$us_assertive_conf_1)==FALSE),],
                   alternative = c("two.sided"), mu=0,
                   paired=FALSE, var.equal=FALSE, conf.level = 0.95) ## t_test is pipe-friendly and the output can be added to a plot using the ggpubr package 

d2$t_polarization<- as.factor(d2$t_polarization)

p2 <- ggplot(d2[which(is.na(d2$us_assertive_conf_1)==FALSE),], aes(x=t_polarization, y=us_assertive_conf_1)) + 
  geom_boxplot(width=0.5,color="grey60", fill="grey", alpha=0.2)+
  stat_summary(fun=mean, geom="point", shape=17,size=5, color="black")+
  geom_jitter(width=0.2,size=1, color="grey60")+
  theme_bw()+
  ylab("Confidence about predictations \n of US assertiveness")+
  xlab("")+
  scale_y_continuous(breaks = seq(0,100,by=10))+
  scale_x_discrete()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 15, b = 0, l = 15), size=30, face="bold"),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0), size=20),
        text = element_text(size=30, face="bold"),
        legend.text=element_text(size=25),
        legend.key.size=unit(0.8, "cm"),
        legend.position="")+
  labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p2

# Figure 4 in Appendix
grid.arrange(p1, p2, 
             ncol = 2, 
             nrow = 1)

ftest1<- var.test(us_active_num ~ t_polarization, d2, 
                  alternative = "two.sided")

ftest2<- var.test(us_assertive_num ~ t_polarization, d2, 
                  alternative = "two.sided")

outcome<- c("US Activeness", "US Assertiveness")
var.control<- c(stat.desc(d2[which(d2$t_polarization=="Control"),]$us_active_num)[[12]],
                stat.desc(d2[which(d2$t_polarization=="Control"),]$us_assertive_num)[[12]])
var.treat<- c(stat.desc(d2[which(d2$t_polarization=="Treatment"),]$us_active_num)[[12]],
              stat.desc(d2[which(d2$t_polarization=="Treatment"),]$us_assertive_num)[[12]])

F.test<- c(paste(paste("F",round(ftest1$statistic, digits=3), sep="="), 
                 paste("p-value",round(ftest1$p.value[[1]], digits=3), sep="="),
                 sep=", "),
           paste(paste("F",round(ftest2$statistic, digits=3), sep="="), 
                 paste("p-value",round(ftest2$p.value[[1]], digits=3), sep="="),
                 sep=", ")
)

dftest<- data.frame(outcome, var.control, var.treat, F.test)

# Table 4 in Appendix

kable(dftest, "latex", booktabs = TRUE)

# Appendix I-H: Exploratory Analysis of Heterogenous Treatment Effects ------------------------------

# High Interest in International Affairs
  # 1 = 30 minutes or more each day on reading about foreign affairs

d2$high_interest = ifelse(d2$intaffairs == "30 minutes to an hour" |
                            d2$intaffairs == "More than an hour", 1, 0)
d2$high_interest<- ifelse(d2$intaffairs=="", NA, d2$high_interest)

m1_highinterest <- lm(china_assertive_num ~ high_interest* polarization_treatment, data = d2)
m2_highinterest <- lm(china_assertive_num ~ high_interest* polarization_treatment  + female + age3554 + 
                  age55plus + han + urban_bi + college_higher, data = d2)

# Access to International Internet
  # 1 = respondent gets news from “International Internet” 


d2$int_internet = ifelse(grepl("International internet", d2$intaffairs_source)==TRUE, 1, 0)
d2$int_internet<- ifelse(d2$intaffairs_source=="", NA, d2$int_internet)

m1_intinternet <- lm(china_assertive_num ~ int_internet* polarization_treatment, data = d2)
m2_intinternet <- lm(china_assertive_num ~ int_internet* polarization_treatment  + female + age3554 + 
                        age55plus + han + urban_bi + college_higher, data = d2)

# Internet Primary News Source
  # 1 = respondent gets news the most from either "Domestic Internet" or "International Internet"

d2$internet_primary = ifelse(grepl("International internet", d2$intaffairs_primary)==TRUE |
                               grepl("Domestic internet", d2$intaffairs_primary)==TRUE, 1, 0)
d2$internet_primary<- ifelse(d2$intaffairs_primary=="", NA, d2$internet_primary)


m1_internetprimary <- lm(china_assertive_num ~ internet_primary* polarization_treatment, data = d2)
m2_internetprimary <- lm(china_assertive_num ~ internet_primary* polarization_treatment  + female + age3554 + 
                       age55plus + han + urban_bi + college_higher, data = d2)

# Any Experience in the U.S.
  # 1 = respondent reports any travel to the U.S.

d2$us_travel_any = ifelse(d2$us_travel=="Yes", 1, 0)
d2$us_travel_any<- ifelse(d2$us_travel=="", NA, d2$us_travel_any)

m1_ustravel <- lm(china_assertive_num ~ us_travel_any* polarization_treatment, data = d2)
m2_ustravel  <- lm(china_assertive_num ~ us_travel_any* polarization_treatment  + female + age3554 + 
                           age55plus + han + urban_bi + college_higher, data = d2)

# Significant Experience in the U.S.
  # 1 = respondent reports spending 1 or more years in the U.S.

d2$us_time_significant = ifelse(d2$us_time=="1-3 years" |
                                  d2$us_time == "More than 3 years", 1, 0)
d2$us_time_significant<- ifelse(d2$us_time=="", NA, d2$us_time_significant)


m1_ussignificant <- lm(china_assertive_num ~ us_time_significant* polarization_treatment, data = d2)
m2_ussignificant  <- lm(china_assertive_num ~ us_time_significant* polarization_treatment  + female + age3554 + 
                     age55plus + han + urban_bi + college_higher, data = d2)

# Communist Party Membership
  # 1 = respondent reports belonging to the Communist Party

d2$ccp_member = ifelse(d2$ccp=="Yes", 1, 0)
d2$ccp_member = ifelse(d2$ccp == "Prefer not to answer", NA, d2$ccp_member)
d2$ccp_member = ifelse(d2$ccp == "", NA, d2$ccp_member)


m1_ccp<- lm(china_assertive_num ~ ccp_member* polarization_treatment, data = d2)
m2_ccp<- lm(china_assertive_num ~ ccp_member* polarization_treatment  + female + age3554 + 
                  age55plus + han + urban_bi + college_higher, data = d2)

# Political Ideology 
  # 5-point Likert from 1 (Very Liberal) to 5 (Very Conservative)

d2$pol_views_scale = ifelse(d2$pol_views=="Liberal", 1, NA)
d2$pol_views_scale = ifelse(d2$pol_views=="Somewhat liberal", 2, d2$pol_views_scale)
d2$pol_views_scale = ifelse(d2$pol_views=="Moderate", 3, d2$pol_views_scale)
d2$pol_views_scale = ifelse(d2$pol_views=="Somewhat conservative", 4, d2$pol_views_scale)
d2$pol_views_scale = ifelse(d2$pol_views=="Very conservative", 5, d2$pol_views_scale)
d2$pol_views_scale = ifelse(d2$pol_views=="Don't know" | d2$pol_views == "Prefer not to  answer", NA, d2$pol_views_scale)
d2$pol_views_scale = as.numeric(d2$pol_views_scale)

m1_polviews<- lm(china_assertive_num ~ pol_views_scale* polarization_treatment, data = d2)
m2_polviews<- lm(china_assertive_num ~ pol_views_scale* polarization_treatment  + female + age3554 + 
              age55plus + han + urban_bi + college_higher, data = d2)

# Hawkishness
  # 1 = respondent reports that China uses military force "too little"

d2$hawkish_bi = ifelse(d2$hawkish == "Too little", 1, 0)
d2$hawkish_bi = ifelse(d2$hawkish == "Prefer not to answer" |
                         d2$hawkish == "Don't know", NA, d2$hawkish_bi)
d2$hawkish_bi = ifelse(d2$hawkish=="", NA, d2$hawkish_bi)

m1_hawkish<- lm(china_assertive_num ~ hawkish_bi* polarization_treatment, data = d2)
m2_hawkish<- lm(china_assertive_num ~ hawkish_bi* polarization_treatment  + female + age3554 + 
          age55plus + han + urban_bi + college_higher, data = d2)

# Appendix Table 5: Exploration of Heterogeneous Effects by Access to Internet and Media

stargazer(m1_highinterest, m2_highinterest, m1_intinternet, m2_intinternet, 
          m1_internetprimary, m2_internetprimary, 
          omit = c("female", "age3554", "age55plus", "han", "urban_bi", "college_higher", "Constant"),
          order = c("polarization_treatment"),
          covariate.labels = c("Polarization Treatment", "High Int. News Interest*Treatment",
                               "Access to Int. Internet*Treatment", "Internet Primary News Source*Treatment",
                               "High Int. News Interest", "Access to International Internet", 
                               "Internet Primary News Source"),
          omit.stat = c("f", "rsq", "ser"),
          add.lines = list(c("Controls", "No", "Yes", "No", "Yes", "No", "Yes")),
          dep.var.caption = "Assertiveness of China",
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)
          
# Appendix Table 6: Exploration of Heterogeneous Effects by Prior U.S. Experience

stargazer(m1_ustravel, m2_ustravel, m1_ussignificant, m2_ussignificant,
          omit = c("female", "age3554", "age55plus", "han", "urban_bi", "college_higher", "Constant"),
          order = c("polarization_treatment"),
          covariate.labels = c("Polarization Treatment", "Any US Travel*Treatment",
                               "Significant US Travel*Treatment", "Any US Travel",
                               "Significant US Travel"),
          omit.stat = c("f", "rsq", "ser"),
          add.lines = list(c("Controls", "No", "Yes", "No", "Yes")),
          dep.var.caption = "Assertiveness of China",
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)

# Appendix Table 7: Exploration of Heterogeneous Effects by Political Affiliation and Ideology

stargazer(m1_ccp, m2_ccp, m1_polviews, m2_polviews, m1_hawkish, m2_hawkish,
          omit = c("female", "age3554", "age55plus", "han", "urban_bi", "college_higher", "Constant"),
          order = c("polarization_treatment"),
          covariate.labels = c("Polarization Treatment", "Communist Party*Treatment",
                               "Political Ideology*Treatment", "Hawkish*Treatment",
                               "Communist Party", "Political Ideology", 
                               "Hawkish"),
          omit.stat = c("f", "rsq", "ser"),
          add.lines = list(c("Controls", "No", "Yes", "No", "Yes", "No", "Yes")),
          dep.var.caption = "Assertiveness of China",
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)

# Appenidx Table 8: Treatment Effects for Hawk and Non-Hawk Respondents

m1_nonhawk<- lm(china_assertive_num ~ polarization_treatment, data = d2[which(d2$hawkish_bi==0),])
summary(m1)

m2_nonhawk<- lm(china_assertive_num ~ polarization_treatment  + female + age3554 + 
          age55plus + han + urban_bi + college_higher, data = d2[which(d2$hawkish_bi==0),])
summary(m2)

m3_hawk<- lm(china_assertive_num ~ polarization_treatment, data = d2[which(d2$hawkish_bi==1),])
summary(m3)

m4_hawk<- lm(china_assertive_num ~ polarization_treatment  + female + age3554 + 
          age55plus + han + urban_bi + college_higher, data = d2[which(d2$hawkish_bi==1),])
summary(m4)

stargazer(m1_nonhawk, m2_nonhawk, m3_hawk, m4_hawk,
          omit = c("female", "age3554", "age55plus", "han", "urban_bi", "college_higher", "Constant"),
          order = c("polarization_treatment"),
          covariate.labels=c("Treatment"),
          omit.stat = c("f", "rsq", "ser"),
          add.lines = list(c("Controls", "No", "Yes", "No", "Yes")),
          dep.var.caption = "Assertiveness of China",
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)

# Manuscript: Observational Analysis --------------------------------------------------
rm(list = ls())

## Import observational data 

dtest<- read.csv("US-China Polarization - Observational Data - Main Sample.csv")

# Figure 4 ----------------------------------------------------------------

## treat only the Rival-US directed dyads as treatment group
dtest$rival<- ifelse(dtest$stateb==2, 1, 0)

## drop dyads without interactions in the month before 2020-11-07
dtest$date<- as.Date(dtest$date, format = "%m/%d/%y")
dpresixty<- filter(dtest, date>=as.Date("2020-10-08", format="%Y-%m-%d") &
                     date<as.Date("2020-11-07", format="%Y-%m-%d"))
dpresixty<- group_by(dpresixty, statea, stateb, id) %>%
  summarise(keep=max(event))

dtest<- left_join(dtest, dpresixty, by=c("statea","stateb","id"))

dtest<- filter(dtest, keep==1)

## break the data into 12 different subsets for different time windows
Jan6<- as.Date("2021-01-06", format="%Y-%m-%d") ## create the reference point (Jan 6th)
dtest$date<- as.Date(dtest$date) ## make sure "date" variable is in Date format
cutpoint<- c(5,10,15,20,25,30,35,40,45,50,55,60) ## create 12 different window width
list_df<- list(length=12) ## create an empty list that will store 12 different data frames (subsamples)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12) 

for (j in number) {
    list_df[[j]]<- filter(dtest, (Jan6-date<=cutpoint[[j]] & Jan6-date>=0) | (date-Jan6<=cutpoint[[j]] & date-Jan6>=0))
} ## use a for-loop to subset the full sample into subsamples with 12 different window width


## DiD Models
list_m<- list(length=12) ## create an empty list that will store 12 different OLS model results

for (j in number) {
  list_m[[j]]<- lm(intensity2 ~ rival*postshock+conttype+logtrade2020+polity2b+polity2a, 
                   data = list_df[[j]])
} ## use a for-loop to run 12 OLS models with different samples stored in list_df

## Table 9 in the Appendix (Table Corresponds to Figure 4 in Manuscript)
stargazer(list_m[[1]], list_m[[2]], list_m[[3]], list_m[[4]], list_m[[5]], list_m[[6]],
          list_m[[7]], list_m[[8]], list_m[[9]], list_m[[10]], list_m[[11]], list_m[[12]],
          align=TRUE, no.space = TRUE,
          column.sep.width = "-22pt",
          covariate.labels=c("US Target", "Post-Jan 6","Contiguity Type", 
                             "Bilateral Trade", "Target Polity", "Initiator Polity",
                             "US Target*Post-Jan 6" ),
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)


## Plot DiD effect
result3 <- data.frame(matrix(nrow = 12, ncol = 3)) ## create an empty data frame that will store the estimated DiD effect and confidence intervals for 12 different windows
colnames(result3) <- c("DiD_effect", "lower","upper")

for (j in number) {
  result3[j, 1]<- list_m[[j]]$coefficients[[8]]
  result3[j, 2]<- confint(list_m[[j]])[8,1]
  result3[j, 3]<- confint(list_m[[j]])[8,2]
} ## use for-loop to calculate the DiD effect (coefficient of the interaction term) and 95% confidence intervals

result3$window<- c(5,10,15,20,25,30,35,40,45,50,55,60)
result3$type<- ifelse(result3$upper<0 | result3$lower>0, 1, 0)
result3$type<- as.factor(result3$type)


p_did <- ggplot(result3,aes(x = window, y = DiD_effect, color=type))+
  geom_point(size=4)+
  scale_color_manual(values=c('gray','black'))+
  geom_errorbar(aes(ymin = lower, ymax=upper),size=1, width=0.5) +
  geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  ylab("DiD Effect: changes in intensity score") +
  xlab("Time Window (number of days before and after Jan 6, 2021)") +
  scale_x_continuous(breaks = round(seq(0, 65, by = 5),1))+
  theme_bw()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 28, b = 0, l = 0), size=25),
        axis.title.x = element_text(margin = margin(t = 30, r = 0, b = 0, l = 0), size=25),
        plot.title = element_text(size=25, face = "bold"))+
  theme(text = element_text(size=30),legend.title = element_blank(),
        strip.text.x = element_text(size = 30),
        legend.background = element_rect(fill="white",size=0.5, linetype="solid", colour ="black"),
        legend.text=element_text(size=30),
        legend.position=c(2.8,0.8),legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(size = 20, face="italic",hjust = 0,vjust = -3),
        plot.margin = margin(2, 2, 2, 2, "cm"))
#+ labs(caption = "Note: Error bars represent 95% confidence intervals. The estimates that are statistically significant 
#         at the 0.05 level are highlighted in black.")

# Figure 4 in Manuscript
p_did


# Figure 5----------------------------------------------------------------------
result4 <- data.frame(matrix(nrow = 12, ncol = 4))
colnames(result4) <- c("meanvalue", "se","lower","upper")


for (j in number) {
  predict1<- mean(predict(list_m[[j]], newdata=mutate(list_df[[j]], rival=1, postshock=1),
                          type="response", na.rm=TRUE))
  
  QI <- data.frame(prdicted.value = predict1)
  
  B <- coef(list_m[[j]])
  V <- vcovCL(list_m[[j]])
  
  set.seed(789)
  sim.coefs <- mvrnorm(1000, mu=B, Sigma=V)
  
  olsm <- list_m[[j]]
  sim.qi <- apply(sim.coefs, 1, FUN=function(y){
    olsm$coefficients <- y
    pred1<- mean(predict(olsm, newdata=mutate(list_df[[j]], rival=1, postshock=1),type="response", na.rm=TRUE))
    
    return(c(pred1))
  })
  sim.qi <- t(sim.qi)
  
  pred.se <- apply(sim.qi, 1, sd)
  pred.lb <- apply(sim.qi, 1, FUN=function(x){
    quantile(x, .025)
  })
  pred.ub <- apply(sim.qi, 1, FUN=function(x){
    quantile(x, .975)
  })
  QI <- mutate(QI, 
               pred.se = pred.se,
               pred.lb = pred.lb,
               pred.ub = pred.ub)
  
  result4[j, 1:4]<- QI
}

result4$window<- c(5,10,15,20,25,30,35,40,45,50,55,60)
result4$type<- ifelse(result4$upper>=0 & result4$lower<=0, 0, 1)
result4$type<- as.factor(result4$type)

p_intensity <- ggplot(result4,aes(x = window, y = meanvalue))+
  geom_point(size=4)+
  #scale_color_manual(values=c('gray','black'))+
  geom_errorbar(aes(ymin = lower, ymax=upper),size=1, width=0.5) +
  geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  ylab("Predicted level of intensity") +
  xlab("Time Window (number of days after Jan 6, 2021)") +
  scale_x_continuous(breaks = round(seq(0, 65, by = 5),1))+
  scale_y_continuous(breaks = round(seq(-0.6, 0, by = 0.3),1))+
  ggtitle("")+
  theme_bw()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 28, b = 0, l = 0), size=25),
        axis.title.x = element_text(margin = margin(t = 30, r = 0, b = 0, l = 0), size=25),
        plot.title = element_text(size=25, face = "bold"))+
  theme(text = element_text(size=30),legend.title = element_blank(),
        strip.text.x = element_text(size = 30),
        legend.background = element_rect(fill="white",size=0.5, linetype="solid", colour ="black"),
        legend.text=element_text(size=30),
        legend.position=c(2.8,0.8),legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(size = 20, hjust = 0,vjust = -3, face="italic"),
        plot.margin = margin(2, 2, 2, 2, "cm"))

# Figure 4
p_intensity

# Appendix Table 10: OLS, US Proteges vs. Other Target

## Code US proteges
#### NATO
dtest$countryb<- countrycode(dtest$stateb, destination = "country.name",
                             origin="cown")

dtest$nato<- ifelse(dtest$countryb=="Albania", 1, 0)
dtest$nato<- ifelse(dtest$countryb=="Belgium", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Bulgaria", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Canada", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Croatia", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Czechia", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Denmark", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Estonia", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="France", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Germany", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Greece", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Hungary", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Iceland", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Italy", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Latvia", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Lithuania", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Luxembourg", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Montenegro", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Netherlands", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="North Macedonia", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Norway", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Poland", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Portugal", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Romania", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Slovakia", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Slovenia", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Spain", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="Turkey", 1, dtest$nato)
dtest$nato<- ifelse(dtest$countryb=="United Kingdom", 1, dtest$nato)

#### Trio Treaty
dtest$trio<- ifelse(dtest$countryb=="Argentina", 1, 0)
dtest$trio<- ifelse(dtest$countryb=="Brazil", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Chile", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Colombia", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Costa Rica", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Dominican Republic", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="EL Salvador", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Guatemala", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Haiti", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Honduras", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Panama", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Paraguay", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Peru", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Bahamas", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Trinidad and Tobago", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Uruguay", 1, dtest$trio)
dtest$trio<- ifelse(dtest$countryb=="Venezuela", 1, dtest$trio)

#### ANZUS
dtest$anzus<- ifelse(dtest$countryb=="Australia", 1, 0)
dtest$anzus<- ifelse(dtest$countryb=="New Zealand", 1, dtest$anzus)

#### Bilateral treaties
dtest$bitreaty<- ifelse(dtest$countryb=="Philippines", 1, 0)
dtest$bitreaty<- ifelse(dtest$countryb=="South Korea", 1, dtest$bitreaty)
dtest$bitreaty<- ifelse(dtest$countryb=="Japan", 1, dtest$bitreaty)
dtest$bitreaty<- ifelse(dtest$countryb=="United Kingdom", 1, dtest$bitreaty)

#### Taiwan and Israel
dtest$othertreaty<- ifelse(dtest$countryb=="Taiwan", 1, 0)
dtest$othertreaty<- ifelse(dtest$countryb=="Israel", 1, dtest$othertreaty)

dtest$rival<- ifelse(dtest$bilat<=0.25 & dtest$nato==1, 1, dtest$rival)
dtest$rival<- ifelse(dtest$bilat<=0.25 & dtest$trio==1, 1, dtest$rival)
dtest$rival<- ifelse(dtest$bilat<=0.25 & dtest$anzus==1, 1, dtest$rival)
dtest$rival<- ifelse(dtest$bilat<=0.25 & dtest$bitreaty==1, 1, dtest$rival)
dtest$rival<- ifelse(dtest$bilat<=0.25 & dtest$othertreaty==1, 1, dtest$rival)

#### Drop the US
dtest<- filter(dtest, stateb!=2)

## break the data into 12 different subsets for different time windows
Jan6<- as.Date("2021-01-06", format="%Y-%m-%d")
dtest$date<- as.Date(dtest$date)
cutpoint<- c(5,10,15,20,25,30,35,40,45,50,55,60)
list_df<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_df[[j]]<- filter(dtest, (Jan6-date<=cutpoint[[j]] & Jan6-date>=0) | (date-Jan6<=cutpoint[[j]] & date-Jan6>=0))
}


## DiD Models
list_m<- list(length=12)

for (j in number) {
  list_m[[j]]<- lm(intensity2 ~ rival*postshock+conttype+logtrade2020+polity2b+polity2a, 
                   data = list_df[[j]])
}

# Appendix Table 10: OLS, US Proteges vs. Other Target

stargazer(list_m[[1]], list_m[[2]], list_m[[3]], list_m[[4]], list_m[[5]], list_m[[6]],
          list_m[[7]], list_m[[8]], list_m[[9]], list_m[[10]], list_m[[11]], list_m[[12]],
          align=TRUE, no.space = TRUE,
          column.sep.width = "-22pt",
          covariate.labels=c("US Protege", "Post-Jan 6","Contiguity Type", 
                             "Bilateral Trade", "Target Polity", "Initiator Polity",
                             "US Protege*Post-Jan 6" ),
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)


result5 <- data.frame(matrix(nrow = 12, ncol = 3))
colnames(result5) <- c("DiD_effect", "lower","upper")

for (j in number) {
  result5[j, 1]<- list_m[[j]]$coefficients[[8]]
  result5[j, 2]<- confint(list_m[[j]])[8,1]
  result5[j, 3]<- confint(list_m[[j]])[8,2]
}

result5$window<- c(5,10,15,20,25,30,35,40,45,50,55,60)
result5$type<- ifelse(result5$upper<0 | result5$lower>0, 1, 0)
result5$type<- as.factor(result5$type)


p_did_protege <- ggplot(result5,aes(x = window, y = DiD_effect, color=type))+
  geom_point(size=4)+
  scale_color_manual(values=c('gray','black'))+
  geom_errorbar(aes(ymin = lower, ymax=upper),size=1, width=0.5) +
  geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  ylab("DiD Effect: changes in intensity score") +
  xlab("Time Window (number of days before and after Jan 6, 2021)") +
  scale_x_continuous(breaks = round(seq(0, 65, by = 5),1))+
  ggtitle("")+
  theme_bw()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 28, b = 0, l = 0), size=25),
        axis.title.x = element_text(margin = margin(t = 30, r = 0, b = 0, l = 0), size=25))+
  theme(text = element_text(size=30),legend.title = element_blank(),
        strip.text.x = element_text(size = 30),
        legend.background = element_rect(fill="white",size=0.5, linetype="solid", colour ="black"),
        legend.text=element_text(size=30),
        legend.position=c(2.8,0.8),legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(size = 20, face="italic",hjust = 0,vjust = -3),
        plot.margin = margin(2, 2, 2, 2, "cm"))+
  labs(caption = "Note: Error bars represent 95% confidence intervals. The estimates that are statistically significant 
         at the 0.05 level are highlighted in black.")
p_did_protege

# Appendix: Observational Analysis ----------------------------------------

## Appendix II-B, Figure 5: Visual check of parallel trends assumption ----------------------------
rm(list = ls())

## Reload Data
dtest<- read.csv("US-China Polarization - Observational Data - Main Sample.csv")

## code the US target as treatment
dtest$rival<- ifelse(dtest$stateb==2, 1, 0)

## break the data into 12 different subsets for different time windows
Jan6<- as.Date("2021-01-06", format="%Y-%m-%d")
dtest$date<- as.Date(dtest$date, format="%m/%d/%y")
cutpoint<- c(5,10,15,20,25,30,35,40,45,50,55,60)
list_df<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_df[[j]]<- filter(dtest, (Jan6-date<=cutpoint[[j]] & Jan6-date>=0) | (date-Jan6<=cutpoint[[j]] & date-Jan6>=0))
}

number<- c(1,2,3,4,5,6,7,8,9,10,11,12)
list_p<- list(length=12)
list_title<- list("5-Day", "10-Day", "15-Day", "20-Day", "25-Day",
                  "30-Day", "35-Day", "40-Day", "45-Day", "50-Day",
                  "55-Day", "60-Day")

for (j in number) {
  dplot<- group_by(list_df[[j]], rival, date) %>%
    summarise(intensity2=mean(intensity, na.rm=TRUE))
  
  dplot$type<- ifelse(dplot$rival==1, "US Target", "Other Target")
  dplot$type<- as.factor(dplot$type)
  dplot$rival<- as.factor(dplot$rival)
  dplot$date<- as.Date(dplot$date)
  
  list_p[[j]] <- ggplot(dplot, aes(x=date, y=intensity2, shape=type, color=type)) +
    geom_point(size=3) + 
    geom_smooth()+
    ggtitle(list_title[[j]])+
    geom_vline(xintercept = as.Date("2021-01-06", format="%Y-%m-%d"),
               color="red") +
    scale_shape_manual(values=c(16, 17))+ 
    scale_color_manual(values=c('blue', 'red'))+
    xlab("")+
    ylab("Intensity Score")+
    scale_x_date(date_labels = "%b-%d")+
    theme(text = element_text(size=20),legend.title = element_blank(),
          plot.title = element_text(size=30),
          legend.background = element_rect(fill="white",size=0.5, linetype="solid", colour ="black"),
          legend.text=element_text(size=35),
          legend.position="right",legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
          plot.caption = element_text(size = 20, hjust = 0,vjust = -3, face="italic"),
          plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm"))
}

ggarrange(list_p[[1]], list_p[[2]], list_p[[3]], list_p[[4]], list_p[[5]], list_p[[6]],
          list_p[[7]], list_p[[8]], list_p[[9]], list_p[[10]], list_p[[11]], list_p[[12]],
          common.legend = TRUE, legend = "right", hjust=-0.2)

## For code that generates Table 9 and 10 in the Appendix, see codes for Figure 4-5

## Appendix II-D, Table 11: Heckman Models (US Target vs. Other Target) -----------------------------------------------------

## Reload Data
rm(list = ls())

dtest<- read.csv("US-China Polarization - Observational Data - Main Sample.csv")


## treat only the Rival-US directed dyads as treatment group
dtest$rival<- ifelse(dtest$stateb==2, 1, 0)

## break the data into 12 different subsets for different time windows
Jan6<- as.Date("2021-01-06", format="%Y-%m-%d")
dtest$date<- as.Date(dtest$date, format="%m/%d/%y")
cutpoint<- c(5,10,15,20,25,30,35,40,45,50,55,60)
list_df<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_df[[j]]<- filter(dtest, (Jan6-date<=cutpoint[[j]] & Jan6-date>=0) | (date-Jan6<=cutpoint[[j]] & date-Jan6>=0))
}

## Heckman models: Rival-US as treatment group
list_m<- list(length=12)

for (j in number) {
  list_m[[j]]<- selection(event ~ rival*postshock+conttype+logtrade2020+polity2b+polity2a,
                          intensity ~ rival*postshock, 
                          data = list_df[[j]], 
                          method="ml")
}



## Table for the selection 
stargazer(list_m[[1]],list_m[[2]],list_m[[3]],list_m[[4]],list_m[[5]],list_m[[6]],
          list_m[[7]],list_m[[8]],list_m[[9]],list_m[[10]],list_m[[11]],list_m[[12]],
          align=TRUE, no.space = TRUE,
          column.sep.width = "-22pt",
          covariate.labels=c("US Target", "Post-Jan 6","Contiguity Type", 
                             "Bilateral Trade", "Target Polity", "Initiator Polity",
                             "US Target*Post-Jan 6"),
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE,
          selection.equation = TRUE)

## Table for the outcome 
stargazer(list_m[[1]],list_m[[2]],list_m[[3]],list_m[[4]],list_m[[5]],list_m[[6]],
          list_m[[7]],list_m[[8]],list_m[[9]],list_m[[10]],list_m[[11]],list_m[[12]],
          align=TRUE, no.space = TRUE,
          column.sep.width = "-22pt",
          covariate.labels=c("US Target", "Post-Jan 6",
                             "US Target*Post-Jan 6"),
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)

## Appendix II-D, Figure 6: Heckman Models  ----------------------------------------------------

result <- data.frame(matrix(nrow = 12, ncol = 4))
colnames(result) <- c("DiD_effect","se","lower","upper")
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

## First plot the DiD effect for the Probit selection part as (Pr(Treated+Post-Jan6)-Pr(Treated+Pre-Jan6))-(Pr(Control+Post-Jan6)-Pr(Control+Pre-Jan6))
## In this algorithm, control variables are held at their observed values for each observation
## Confidence intervals are gained through simulations of 1000 draws of new coefficients (results will differ slightly from appendix)

for (j in number) {
  simdf<- list_df[[j]]
  simdf<- filter(simdf, is.na(polity2b)==FALSE)
  predict1<- mean(predict(list_m[[j]], newdata=mutate(simdf, rival=1, postshock=1),
                          type="response", part="selection",na.rm=TRUE))
  
  predict2<- mean(predict(list_m[[j]], newdata=mutate(simdf, rival=1, postshock=0),
                          type="response", part="selection",na.rm=TRUE))
  
  predict3<- mean(predict(list_m[[j]], newdata=mutate(simdf, rival=0, postshock=1),
                          type="response", part="selection",na.rm=TRUE))
  
  predict4<- mean(predict(list_m[[j]], newdata=mutate(simdf, rival=0, postshock=0),
                          type="response", part="selection",na.rm=TRUE))
  
  DiD<- (predict1-predict2)-(predict3-predict4)
  
  QI <- data.frame(prdicted.value = DiD)
  
  B <- coef(list_m[[j]])
  V <- vcovCL(list_m[[j]])
  
  sim.coefs <- mvrnorm(1000, mu=B, Sigma=V)
  
  m_sim <- list_m[[j]]
  sim.qi <- apply(sim.coefs, 1, FUN=function(y){
    m_sim$estimate <- y
    predict1<- mean(predict(m_sim, newdata=mutate(simdf, rival=1, postshock=1),
                            type="response", part="selection",na.rm=TRUE))
    
    predict2<- mean(predict(m_sim, newdata=mutate(simdf, rival=1, postshock=0),
                            type="response", part="selection",na.rm=TRUE))
    
    predict3<- mean(predict(m_sim, newdata=mutate(simdf, rival=0, postshock=1),
                            type="response", part="selection",na.rm=TRUE))
    
    predict4<- mean(predict(m_sim, newdata=mutate(simdf, rival=0, postshock=0),
                            type="response", part="selection",na.rm=TRUE))
    
    DiD<- (predict1-predict2)-(predict3-predict4)
    
    return(c(DiD))
  })
  sim.qi <- t(sim.qi)
  
  pred.se <- apply(sim.qi, 1, sd)
  pred.lb <- apply(sim.qi, 1, FUN=function(x){
    quantile(x, .025)
  })
  pred.ub <- apply(sim.qi, 1, FUN=function(x){
    quantile(x, .975)
  })
  QI <- mutate(QI, 
               pred.se = pred.se,
               pred.lb = pred.lb,
               pred.ub = pred.ub)
  
  result[j, 1:4]<- QI
}

result$window<- c(5,10,15,20,25,30,35,40,45,50,55,60)
result$type<- ifelse(result$upper>=0 & result$lower<=0, 0, 1)
result$type<- as.factor(result$type)


g1 <- ggplot(result,aes(x = window, y = DiD_effect, color=type))+
  geom_point(size=4)+
  scale_color_manual(values=c('gray','black'))+
  geom_errorbar(aes(ymin = lower, ymax=upper),size=1, width=0.5) +
  geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  ylab("DiD Effect: changes in the probability of taking any action") +
  xlab("Time Window (number of days before and after Jan 6, 2021)") +
  scale_x_continuous(breaks = round(seq(0, 65, by = 5),1))+
  ggtitle("(a): Changes in the probability of taking any action")+
  theme_bw()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 28, b = 0, l = 0), size=25, face="plain"),
        axis.title.x = element_text(margin = margin(t = 30, r = 0, b = 0, l = 0), size=25, face="plain"),
        title=element_text(size=25, face="bold"))+
  theme(text = element_text(size=30),legend.title = element_blank(),
        strip.text.x = element_text(size = 30),
        legend.background = element_rect(fill="white",size=0.5, linetype="solid", colour ="black"),
        legend.text=element_text(size=30),
        legend.position=c(2.8,0.8),legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        #plot.caption = element_text(size = 18, hjust = 0,vjust = -3, face="italic", color = "white"),
        plot.margin = margin(2, 2, 2, 2, "cm"))
g1

## Plot the DiD effect on outcome
result2 <- data.frame(matrix(nrow = 12, ncol = 4))
colnames(result2) <- c("DiD_effect", "se", "lower","upper")
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

## Second plot the DiD effect for the outcome part using the same algorithm
for (j in number) {
  simdf<- list_df[[j]]
  simdf<- filter(simdf, is.na(polity2b)==FALSE)
  predict1<- mean(predict(list_m[[j]], newdata=mutate(simdf, rival=1, postshock=1),
                          type="conditional", part="outcome",na.rm=TRUE))
  
  predict2<- mean(predict(list_m[[j]], newdata=mutate(simdf, rival=1, postshock=0),
                          type="conditional", part="outcome",na.rm=TRUE))
  
  predict3<- mean(predict(list_m[[j]], newdata=mutate(simdf, rival=0, postshock=1),
                          type="conditional", part="outcome",na.rm=TRUE))
  
  predict4<- mean(predict(list_m[[j]], newdata=mutate(simdf, rival=0, postshock=0),
                          type="conditional", part="outcome",na.rm=TRUE))
  
  DiD<- (predict1-predict2)-(predict3-predict4)
  
  QI <- data.frame(prdicted.value = DiD)
  
  B <- coef(list_m[[j]])
  V <- vcovCL(list_m[[j]])
  
  sim.coefs <- mvrnorm(1000, mu=B, Sigma=V)
  
  m_sim <- list_m[[j]]
  sim.qi <- apply(sim.coefs, 1, FUN=function(y){
    m_sim$estimate <- y
    predict1<- mean(predict(m_sim, newdata=mutate(simdf, rival=1, postshock=1),
                            type="conditional", part="outcome",na.rm=TRUE))
    
    predict2<- mean(predict(m_sim, newdata=mutate(simdf, rival=1, postshock=0),
                            type="conditional", part="outcome",na.rm=TRUE))
    
    predict3<- mean(predict(m_sim, newdata=mutate(simdf, rival=0, postshock=1),
                            type="conditional", part="outcome",na.rm=TRUE))
    
    predict4<- mean(predict(m_sim, newdata=mutate(simdf, rival=0, postshock=0),
                            type="conditional", part="outcome",na.rm=TRUE))
    
    DiD<- (predict1-predict2)-(predict3-predict4)
    
    return(c(DiD))
  })
  sim.qi <- t(sim.qi)
  
  pred.se <- apply(sim.qi, 1, sd)
  pred.lb <- apply(sim.qi, 1, FUN=function(x){
    quantile(x, .025)
  })
  pred.ub <- apply(sim.qi, 1, FUN=function(x){
    quantile(x, .975)
  })
  QI <- mutate(QI, 
               pred.se = pred.se,
               pred.lb = pred.lb,
               pred.ub = pred.ub)
  
  result2[j, 1:4]<- QI
}

result2$window<- c(5,10,15,20,25,30,35,40,45,50,55,60)
result2$type<- ifelse(result2$upper<0 | result2$lower>0, 1, 0)
result2$type<- as.factor(result2$type)


g2 <- ggplot(result2,aes(x = window, y = DiD_effect, color=type))+
  geom_point(size=4)+
  scale_color_manual(values=c('gray','black'))+
  geom_errorbar(aes(ymin = lower, ymax=upper),size=1, width=0.5) +
  geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  ylab("DiD Effect: changes in intensity score") +
  xlab("Time Window (number of days before and after Jan 6, 2021)") +
  scale_x_continuous(breaks = round(seq(0, 65, by = 5),1))+
  ggtitle("(b): Changes in intensity score")+
  theme_bw()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 28, b = 0, l = 0), size=25, face="plain"),
        axis.title.x = element_text(margin = margin(t = 30, r = 0, b = 0, l = 0), size=25, face="plain"),
        title=element_text(size=25, face="bold"))+
  theme(text = element_text(size=30),legend.title = element_blank(),
        strip.text.x = element_text(size = 30),
        legend.background = element_rect(fill="white",size=0.5, linetype="solid", colour ="black"),
        legend.text=element_text(size=30),
        legend.position=c(2.8,0.8),legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        #plot.caption = element_text(size = 20, hjust = 0,vjust = -3, face="italic", color="white"),
        plot.margin = margin(2, 2, 2, 2, "cm"))
g2

g3<- ggarrange(g1, g2, ncol = 2, nrow = 1)
g3+theme(plot.caption = element_text(size = 20, hjust = 0.1,vjust = -3),
         plot.margin = margin(2, 2, 2, 2, "cm"))+
  labs(caption = "Note: Error bars represent 95% confidence intervals generated through 1000 draws of new coefficients. \n    The estimates that are statistically significant at the 0.05 level are highlighted in black")

## Appendix II-E, Table 12: OLS Models with SEs Clustered on Directed-Dyads: U.S. Target vs. Other Target --------------------------------

## Reload data
rm(list = ls())

dtest<- read.csv("US-China Polarization - Observational Data - Main Sample.csv")

## drop dyads without interactions in the month before 2020-11-07
dtest$date<- as.Date(dtest$date, format="%m/%d/%y")
dpresixty<- filter(dtest, date>=as.Date("2020-10-08", format="%Y-%m-%d") &
                     date<as.Date("2020-11-07", format="%Y-%m-%d"))
dpresixty<- group_by(dpresixty, statea, stateb, id) %>%
  summarise(event=max(event))

dpresixty<- dplyr::rename(dpresixty, keep=event)

dtest<- left_join(dtest, dpresixty, by=c("statea","stateb","id"))

dtest<- filter(dtest, keep==1)

## treat only the Rival-US directed dyads as treatment group
dtest$rival<- ifelse(dtest$stateb==2, 1, 0)

## break the data into 12 different subsets for different time windows
Jan6<- as.Date("2021-01-06", format="%Y-%m-%d")
cutpoint<- c(5,10,15,20,25,30,35,40,45,50,55,60)
list_df<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_df[[j]]<- filter(dtest, (Jan6-date<=cutpoint[[j]] & Jan6-date>=0) | (date-Jan6<=cutpoint[[j]] & date-Jan6>=0))
}


list_m<- list(length=12)
list_se<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_m[[j]]<- lm(intensity2 ~ rival*postshock+conttype+logtrade2020+polity2b+polity2a, 
                   data = list_df[[j]])
  list_se[[j]]<- coeftest(list_m[[j]],vcov = vcovCL(list_m[[j]], cluster=list_df[[j]]$dyad))[,2]
}

stargazer(list_m[[1]], list_m[[2]], list_m[[3]], list_m[[4]], list_m[[5]], list_m[[6]],
          list_m[[7]], list_m[[8]], list_m[[9]], list_m[[10]], list_m[[11]], list_m[[12]],
          se=list(list_se[[1]],list_se[[2]],list_se[[3]],list_se[[4]],list_se[[5]],list_se[[6]],
                  list_se[[7]],list_se[[8]],list_se[[9]],list_se[[10]],list_se[[11]],list_se[[12]]),
          align=TRUE, no.space = TRUE,
          column.sep.width = "-22pt",
          covariate.labels=c("US Target", "Post-Jan 6","Contiguity Type", 
                             "Bilateral Trade", "Target Polity", "Initiator Polity",
                             "US Target*Post-Jan 6" ),
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)


## Appendix II-E, Table 13: OLS, Two-Period DiD, U.S. Target vs. Other Target --------------------------------

## Reload Data
rm(list = ls())

dtest<- read.csv("US-China Polarization - Observational Data - Main Sample.csv")

## treat only the Rival-US directed dyads as treatment group
dtest$rival<- ifelse(dtest$stateb==2, 1, 0)

## Drop directed dyads without interactions in the month before 2020-11-07
dtest$date<- as.Date(dtest$date, format="%m/%d/%y")
dpresixty<- filter(dtest, date>=as.Date("2020-10-08", format="%Y-%m-%d") &
                     date<as.Date("2020-11-07", format="%Y-%m-%d"))
dpresixty<- group_by(dpresixty, statea, stateb, id) %>%
  summarise(event=max(event))

dpresixty<- dplyr::rename(dpresixty, keep=event)

dtest<- left_join(dtest, dpresixty, by=c("statea","stateb","id"))

dtest<- filter(dtest, keep==1)

dtest$treat<- ifelse(dtest$postshock==1 & dtest$rival==1, 1, 0)
dtest$dyad<- paste(dtest$statea, dtest$stateb, sep = "-")

## break the data into 12 different subsets for different time windows
## collapse each subsample into pre- and post-Jan 6 two-period
Jan6<- as.Date("2021-01-06", format="%Y-%m-%d")
cutpoint<- c(5,10,15,20,25,30,35,40,45,50,55,60)
list_df<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_df[[j]]<- filter(dtest, (Jan6-date<=cutpoint[[j]] & Jan6-date>=0) | (date-Jan6<=cutpoint[[j]] & date-Jan6>=0))
  list_df[[j]]<- group_by(list_df[[j]], statea, stateb, postshock, treat,  dyad,  rival) %>%
    summarize(intensity2=mean(intensity2, na.rm=TRUE), 
              conttype=mean(conttype, na.rm=TRUE),
              logtrade2020=mean(logtrade2020, na.rm=TRUE),
              polity2b=mean(polity2b, na.rm=TRUE),
              polity2a=mean(polity2a, na.rm=TRUE))
}


list_m<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_m[[j]]<- lm(intensity2 ~ rival*postshock+conttype+logtrade2020+polity2b+polity2a, 
                   data = list_df[[j]])
}

stargazer(list_m[[1]], list_m[[2]], list_m[[3]], list_m[[4]], list_m[[5]], list_m[[6]],
          list_m[[7]], list_m[[8]], list_m[[9]], list_m[[10]], list_m[[11]], list_m[[12]],
          align=TRUE, no.space = TRUE,
          column.sep.width = "-22pt",
          covariate.labels=c("US Target", "Post-Jan 6","Contiguity Type", 
                             "Bilateral Trade", "Target Polity", "Initiator Polity",
                             "US Target*Post-Jan 6" ),
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)

## Appendix II-E, Table 14: OLS, Only Hostile Interactions, U.S. Target vs. Other Target --------------------------------

## Reload Data
rm(list = ls())

dtest<- read.csv("US-China Polarization - Observational Data - Main Sample.csv")

## treat only the Rival-US directed dyads as treatment group
dtest$rival<- ifelse(dtest$stateb==2, 1, 0)

## mark directed dyads without interactions in the month before 2020-11-07
dtest$date<- as.Date(dtest$date, format="%m/%d/%y")
dpresixty<- filter(dtest, date>=as.Date("2020-10-08", format="%Y-%m-%d") &
                     date<as.Date("2020-11-07", format="%Y-%m-%d"))
dpresixty<- group_by(dpresixty, statea, stateb, id) %>%
  summarise(event=max(event))

dpresixty<- dplyr::rename(dpresixty, keep=event)

dtest<- left_join(dtest, dpresixty, by=c("statea","stateb","id"))

dtest<- filter(dtest, keep==1)


## break the data into 12 different subsets for different time windows
Jan6<- as.Date("2021-01-06", format="%Y-%m-%d")
dtest$date<- as.Date(dtest$date)
cutpoint<- c(5,10,15,20,25,30,35,40,45,50,55,60)
list_df<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_df[[j]]<- filter(dtest, (Jan6-date<=cutpoint[[j]] & Jan6-date>=0) | (date-Jan6<=cutpoint[[j]] & date-Jan6>=0))
}




list_m<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_m[[j]]<- lm(negintense2 ~ rival*postshock+conttype+logtrade2020+polity2b+polity2a, 
                   data = list_df[[j]])
}

stargazer(list_m[[1]], list_m[[2]], list_m[[3]], list_m[[4]], list_m[[5]], list_m[[6]],
          list_m[[7]], list_m[[8]], list_m[[9]], list_m[[10]], list_m[[11]], list_m[[12]],
          align=TRUE, no.space = TRUE,
          column.sep.width = "-22pt",
          covariate.labels=c("US Target", "Post-Jan 6","Contiguity Type", 
                             "Bilateral Trade", "Target Polity", "Initiator Polity",
                             "US Target*Post-Jan 6" ),
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)

## Appendx II-E, Table 15: OLS, Democratic Rival vs. Other Target (Spillover) ---------------------
## Reload Data
rm(list = ls())

dtest<- read.csv("US-China Polarization - Observational Data - Main Sample.csv")

## treat only US Rivals' own democratic rivals as treatment group (to examine possible spillovers)
dtest$rival<- ifelse(dtest$bilat<=0.25 & dtest$polity2b>=6, 1, 0)

## drop the US
dtest<- filter(dtest, stateb!=2)

## mark directed dyads without interactions in the month before 2020-11-07
dtest$date<- as.Date(dtest$date, format="%m/%d/%y")
dpresixty<- filter(dtest, date>=as.Date("2020-10-08", format="%Y-%m-%d") &
                     date<as.Date("2020-11-07", format="%Y-%m-%d"))
dpresixty<- group_by(dpresixty, statea, stateb, id) %>%
  summarise(event=max(event))

dpresixty<- dplyr::rename(dpresixty, keep=event)

dtest<- left_join(dtest, dpresixty, by=c("statea","stateb","id"))

dtest<- filter(dtest, keep==1)

## break the data into 12 different subsets for different time windows
Jan6<- as.Date("2021-01-06", format="%Y-%m-%d")
cutpoint<- c(5,10,15,20,25,30,35,40,45,50,55,60)
list_df<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_df[[j]]<- filter(dtest, (Jan6-date<=cutpoint[[j]] & Jan6-date>=0) | (date-Jan6<=cutpoint[[j]] & date-Jan6>=0))
}


list_m<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_m[[j]]<- lm(intensity2 ~ rival*postshock+conttype+logtrade2020+polity2b+polity2a, 
                   data = list_df[[j]])
}

stargazer(list_m[[1]], list_m[[2]], list_m[[3]], list_m[[4]], list_m[[5]], list_m[[6]],
          list_m[[7]], list_m[[8]], list_m[[9]], list_m[[10]], list_m[[11]], list_m[[12]],
          align=TRUE, no.space = TRUE,
          column.sep.width = "-22pt",
          covariate.labels=c("Democratic Rivals", "Post-Jan 6","Contiguity Type", 
                             "Bilateral Trade", "Target Polity", "Initiator Polity",
                             "Democratic Rivals*Post-Jan 6" ),
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)

## Appendix II-E, Table 16: OLS, U.S. Target vs. Other Target, Politically Relevant Dyad Sample -----------------
## Reload Data
rm(list = ls())

dtest<- read.csv("US-China Polarization - Observational Data - PRD Sample.csv")

## Treat Rival-US dyads as treated
dtest$rival<- ifelse(dtest$stateb==2, 1, 0)


## drop dyads without interactions in the month before 2020-11-07
dtest$date<- as.Date(dtest$date, format="%m/%d/%y")
dpresixty<- filter(dtest, date>=as.Date("2020-10-08", format="%Y-%m-%d") &
                     date<as.Date("2020-11-07", format="%Y-%m-%d"))
dpresixty<- group_by(dpresixty, statea, stateb) %>%
  summarise(event=max(event))

dpresixty<- dplyr::rename(dpresixty, keep=event)

dtest<- left_join(dtest, dpresixty, by=c("statea","stateb"))

dtest<- filter(dtest, keep==1)

## break the data into 12 different subsets for different time windows
Jan6<- as.Date("2021-01-06", format="%Y-%m-%d")
cutpoint<- c(5,10,15,20,25,30,35,40,45,50,55,60)
list_df<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_df[[j]]<- filter(dtest, (Jan6-date<=cutpoint[[j]] & Jan6-date>=0) | (date-Jan6<=cutpoint[[j]] & date-Jan6>=0))
}


## DiD Models
list_m<- list(length=12)

for (j in number) {
  list_m[[j]]<- lm(intensity2 ~ rival*postshock+conttype+logtrade2020+polity2b+polity2a, 
                   data = list_df[[j]])
}

stargazer(list_m[[1]], list_m[[2]], list_m[[3]], list_m[[4]], list_m[[5]], list_m[[6]],
          list_m[[7]], list_m[[8]], list_m[[9]], list_m[[10]], list_m[[11]], list_m[[12]],
          align=TRUE, no.space = TRUE,
          column.sep.width = "-22pt",
          covariate.labels=c("US Target", "Post-Jan 6","Contiguity Type", 
                             "Bilateral Trade", "Target Polity", "Initiator Polity",
                             "US Target*Post-Jan 6" ),
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)

## Appendix II-E, Table 17: OLS, U.S. Rivals–U.S. vs. Non-Rivals–U.S (Alternative Comparison Group) -----------------

## Reload Data
rm(list = ls())

dtest<- read.csv("US-China Polarization - Observational Data - Only US Target Sample.csv")

## Code U.S. Rival-U.S. as treated
dtest$rival<- ifelse(dtest$rival_pd==1, 1, 0)

## drop dyads without interactions in the month before 2020-11-07
dtest$date<- as.Date(dtest$date, format="%m/%d/%y")
dpresixty<- filter(dtest, date>=as.Date("2020-10-08", format="%Y-%m-%d") &
                     date<as.Date("2020-11-07", format="%Y-%m-%d"))
dpresixty<- group_by(dpresixty, ccode) %>%
  summarise(event=max(event))

dpresixty<- dplyr::rename(dpresixty, keep=event)

dtest<- left_join(dtest, dpresixty, by=c("ccode"))

dtest<- filter(dtest, keep==1)

## break the data into 12 different subsets for different time windows
Jan6<- as.Date("2021-01-06", format="%Y-%m-%d")
cutpoint<- c(5,10,15,20,25,30,35,40,45,50,55,60)
list_df<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_df[[j]]<- filter(dtest, (Jan6-date<=cutpoint[[j]] & Jan6-date>=0) | (date-Jan6<=cutpoint[[j]] & date-Jan6>=0))
}


## DiD Models
list_m<- list(length=12)

for (j in number) {
  list_m[[j]]<- lm(intensity2 ~ rival*postshock+conttype+logtrade2020+polity2, 
                   data = list_df[[j]])
}

## Appendix Table 17 
stargazer(list_m[[1]], list_m[[2]], list_m[[3]], list_m[[4]], list_m[[5]], list_m[[6]],
          list_m[[7]], list_m[[8]], list_m[[9]], list_m[[10]], list_m[[11]], list_m[[12]],
          align=TRUE, no.space = TRUE,
          column.sep.width = "-22pt",
          covariate.labels=c("US Rival Initiator", "Post-Jan 6","Contiguity Type", 
                             "Bilateral Trade", "Initiator Polity",
                             "US Target*Post-Jan 6" ),
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)

## Appendix II-E, Table 18: OLS, NATO Members-U.S. vs. NATO Members-Others (Allies Sample) -----------------

## Reload Data
rm(list = ls())

dtest<- read.csv("US-China Polarization - Observational Data - NATO Initiator Sample.csv")

## drop dyads without interactions in the month before 2020-11-07
dtest$date<- as.Date(dtest$date, format="%m/%d/%y")
dpresixty<- filter(dtest, date>=as.Date("2020-10-08", format="%Y-%m-%d") &
                     date<as.Date("2020-11-07", format="%Y-%m-%d"))
dpresixty<- group_by(dpresixty, statea, stateb) %>%
  summarise(event=max(event))

dpresixty<- dplyr::rename(dpresixty, keep=event)

dtest<- left_join(dtest, dpresixty, by=c("statea","stateb"))

dtest<- filter(dtest, keep==1)

## break the data into 12 different subsets for different time windows
Jan6<- as.Date("2021-01-06", format="%Y-%m-%d")
cutpoint<- c(5,10,15,20,25,30,35,40,45,50,55,60)
list_df<- list(length=12)
number<- c(1,2,3,4,5,6,7,8,9,10,11,12)

for (j in number) {
  list_df[[j]]<- filter(dtest, (Jan6-date<=cutpoint[[j]] & Jan6-date>=0) | (date-Jan6<=cutpoint[[j]] & date-Jan6>=0))
}


## DiD Models
list_m<- list(length=12)

for (j in number) {
  list_m[[j]]<- lm(intensity2 ~ treat*postshock+conttype+logtrade2020+polity2b+polity2a, 
                   data = list_df[[j]])
}

## Appendix Table 18
stargazer(list_m[[1]], list_m[[2]], list_m[[3]], list_m[[4]], list_m[[5]], list_m[[6]],
          list_m[[7]], list_m[[8]], list_m[[9]], list_m[[10]], list_m[[11]], list_m[[12]],
          align=TRUE, no.space = TRUE,
          column.sep.width = "-22pt",
          covariate.labels=c("US Target", "Post-Jan 6","Contiguity Type", 
                             "Bilateral Trade", "Target Polity", "Initiator Polity",
                             "US Target*Post-Jan 6" ),
          star.char = c("+", "*", "**", "***"), 
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.align = "r",
          notes.append = FALSE,
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001",
          dep.var.labels.include = FALSE)

